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Abstract 

In this paper, 14-velocity and 18-velocity multiple-relaxation-time (MRT) 
lattice Boltzmann (LB) models are proposed for three-dimensional incom¬ 
pressible flows. These two models are constructed based on the incompress¬ 
ible LBGK model proposed by He et ah (Chin. Phys., 2004, 13: 40-46) and 
the MRT LB model proposed by d’Humieres et ah (Philos. Trans. R. Soc., 
A, 2002, 360: 437-451), which have advantages in the computational effi¬ 
ciency and stability, respectively. Through the Chapman-Enskog analysis, 
the models can recover to three-dimensional incompressible Navier-Stokes 
equations in the low Mach number limit. To verify the present models, the 
steady Poiseuille flow, unsteady pulsatile flow and lid-driven cavity flow in 
three dimensions are simulated. The simulation results agree well with the 
analytical solutions or the existing numerical results. Moreover, it is found 
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that the present models show higher accuracy than d’Humieres et ah model 
and better stability than He et ah model. 

Keywords: lattice Boltzmann model, multiple-relaxation-time, 
three-dimensional incompressible flow, pulsatile flow in 3D rectangular 
channel 


1. Introduction 


The lattice Boltzmann method is an innovative a pp roach based on ki¬ 
netic theory to simulate various complex fluid systems [l|, l2| . The signiflcant 
advantages of lattice Boltzmann method are the natural parallelism of algo¬ 
rithm, simplicity of programming and ease of dealing with complex boundary 
conditions. It has been successfully applied in the held of complex fluids, such 
a. .nuMphaae EuidaS. ...oEuidica g fluids u. pouous n.edia gg and 
impinging fluids P, [8 . 

Until now, the lattice Bhatnagar-Gross-Krook (LBGK) model, which is 
based on a single-relaxation-time (SRT) approximation jo], is still the most 
popular LB model due to its extreme simplicity. The earliest LBGK model is 
proposed by Qian et al. [1^, which is often used to simulate the incompress¬ 
ible flow in the low Mach number limit. However, through the Ghapman- 
Enskog (G-E) procedure, this model can only recover to the compressible 
Navier-Stokes (N-S) equations in the low Mach number limit. If the density 
fluctuation is assumed to be negligible, the incompressible N-S equations can 
be derived. But in practical simulations, sometimes the density fluctuation 
cannot be neglected. In this case, there is compressible effect in the simula¬ 
tions and this effect may lead to serious numerical errors. In fact, Qian et al. 
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model can be viewed as a compressible scheme to simulate incompressible 
flows. There are efforts to weaken the low Mach number restriction of Qian 
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, ll2| , while in our 


et ah model to extend this model for compressible flows 
paper we are focused on how to reduce the compressible effect in existing LB 
models. 

In order to reduce the compressible effect in Qian et al. model, many other 
LB models have been proposed 13l4l6l|. Among them, the models proposed 
by He and Luo, and Guo et al. are widely used. The basic idea of He-Luo 


model 


14| is to neglect the terms of higher order Mach number in equilibrium 


distribution function, which can explicitly reduce the compressible effect as 
demonstrated in the following simulations. However, He-Luo model can only 
accurately recover to the incompressible momentum equations, but keeps the 
term -^dp/dt in the continuity equation, where p = (?p /po is the normalized 
pressure. When He-Luo model is applied to the unsteady flow, it requires an 
additional condition T S> L/cs (T and L are characteristic time and length, 
respectively), to eliminate the compressible effect. 

As we know, the incompressible limit is equivalent to low Mach number 
limit. To overcome the shortcoming of He-Luo model, Guo et al. proposed a 
new LBGK model for two-dimensional incompressible flows. Guo et al. 
model can exactly recover to the incompressible N-S equations only in the 
low Mach number limit, which is accomplished by completely decoupling the 
pressure and density and delicately designing the equilibrium distribution 
function. To our knowledge, Guo et al. model is the first LBGK model 
which can be applied to steady and unsteady incompressible flows while 
simultaneously eliminating the compressible effect completely. Due to the 
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advantage of Guo et al. model, He et al. extended this model to three 
dimensions and proposed the three-dimensional 15-velocity and 19-velocity 


LBGK models for incompressible flows 


16|. 


Although the LBGK model has a very simple algorithm and is popularly 
used, its stability is not always satisfactory in the practical applications. To 
overcome this shortcoming, many other LB models have been developed in 


the past few years 17H24|. Among them, the MRT LB model 


3 - 24 | 


has re¬ 


ceived the most attention due to its superior numerical stability. d’Humieres 


hrstly proposed the MRT LB model 


221 ] at nearly the same time with Qian 


et al. mode 
of model 


Lallemand and Luo carried out detailed analysis on this type 


23| and found that it has much better performance than the LBGK 


model in the stability. To further demonstrate the superior stability of MRT 
model over the LBGK model, d’Humieres et al. developed three-dimensional 
15-velocity and 19-velocity MRT models 241. 

The MRT model has much better stability than the LBGK model, but 
in the aspect of computational efficiency the MRT model could be about 
15% slower than the LBGK counterpart in terms of the number of nodes 
updated per second [2^. The update of one node includes the memory 


access and the floating-point operations, so the computational efficiency of 
MRT LB schemes is mostly limited by the memory access quantity and the 
calculation amount. To improve the computational efficiency of the present 
MRT models, we propose a new class of three-dimensional MRT models with 
less memory consumption and calculation amount in this paper. 

The MRT model proposed by d’Humieres et al. for three-dimensional 


incompressible flows 


24| is based on the Qian et al. model or He-Luo model. 
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which both use q discrete velocity directions in d dimensions, i.e., use the 
DdQg lattice models. However, it is noticed that the incompressible LBGK 


model proposed by He et ah 


16| only uses q — 1 discrete velocity directions 


m 


the actual computation. Moreover, d’Humieres et al. MRT model contains 
a moment corresponding to kinetic energy square, which does not affect the 
recovered macroscopic N-S equations. Therefore, based on the He et al. 
model and d’Humieres et al. model, it is possible to construct an MRT 
model with a (g — 1) x (g — 1) transformation matrix, which can reduce 
the memory consumption and enhance computation efficiency of the existing 
MRT models in three dimensions. 

The above idea is enlightened by the work of Du and Shi, who proposed a 
two-dimensional 8-velocity incompressible (iD2Q8) MRT model 2^ based on 
Guo et al. LBGK model and the two-dimensional MRT model proposed by 
Lallemand and Luo. As a continuing work, we propose two three-dimensional 
MRT models with 14-velocity and 18-velocity based on He et al. LBGK 
model and d’Humieres et al. MRT model in this paper. The general con¬ 
struction method of (g — 1) x (g — 1) transformation matrix is presented. 
Through the G-E procedure, the proposed models can recover to the incom¬ 
pressible N-S equations in the low Mach number limit. The numerical results 
of unsteady pulsatile flow and cavity flow show that the proposed model is 
more accurate than d’Humieres et al. MRT model and more stable than He 
et al. LBGK model, where d’Humieres et al. MRT model and He et al. 
LBGK model are two widely used LB models for three-dimensional incom¬ 
pressible flows. As an example, only the 14-velocity model is presented in 
details in this paper. 
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The rest of the paper is organized as follows. Section [2] briefly describes 
the three-dimensional 15-velocity incompressible (iD3Q15) LBGK model pro¬ 
posed by He et ah. Section |3] presents onr three-dimensional 14-velocity in¬ 
compressible (iD3Q14) MRT model. We provides the simnlation resnlts for 
three benchmark problems: the three-dimensional Poisenille flow, pnlsatile 
flow and cavity flow by nsing the proposed iD3Q14 MRT model in section 
m We also compared some resnlts with those coming from d’Hnmieres et 
al. D3Q15 MRT model and He et ah iD3Q15 LBGK model. Section [5] con- 


clndes this paper. Appendix A briefly give the derivation of incompressible 
N-S eqnations from iD3Q14 MRT model. Appendix B ontlines the three- 
dimensional 18-velocity incompressible (iD3Q18) MRT model. 


2. iD3Q15 LBGK model proposed by He et al. 


The iD3Q15 LBGK model proposed by He et al. in Ref. 
discrete velocity directions as follows: 



inclndes 15 


{Co, Cl,. . . , C14} 


^ 0 1 -1 0 0 0 0 1 -1 1 -1 1 -1 

0001 -1 001 1 -1 -1 1 1 

^ 000001-111 1 1 - 1-1 


1 -1 
-1 -1 
-1 -1 


where c = is the particle velocity and bx and bt are the lattice spacing 
and time step, respectively. 

The evolntion eqnation of the dynamical system is 


/„(a; ++ - f„{x,t) = --{fo,{x,t) - f^^\x,t)) (1) 

T 

where fa{x^t) and fa'^\x,t) are the distribntion fnnction and eqnilibrinm 
distribntion fnnction of particle with velocity Cq at node x and time t, r is 
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dimensionless relaxation time. For the iD3Q15 LBGK model, the equilibrium 
distribution function is designed as 


J a 


Aq + UJa{ -5-h 


(c„ ■ uf 

2cf 


\u\ 


■), 


( 2 ) 


where 


and 


A. 


Po + (wo - l)p/c^, a = 0, 

Uap/cl, a = 1 — 18, 




1/3, 

< 1/18, 
^ 1/36, 


a = 0, 
a = 1 — 6, 
a = 7 - 18, 


(3) 


(4) 


Cs = cj \/3 is the sound speed, p and u are the pressure and velocity, po is a 
constant. 


The macroscopic flow velocity u and pressure p are obtained by: 


14 

It = ^ Cafa, (5a) 

a=l 


(5b) 

a=l * 

From Eqs. (j5]), we can see that the computation of macroscopic quantities 
only require the distribution functions in 14 discrete velocity directions, so 
iD3Q15 LBGK model is actually a 14-velocity incompressible LBGK model. 

Through the G-E expansion, the incompressible Navier-Stokes equations 
can be recovered as 

V • tt = 0, 

+ u ■ Vu = —Vp -|- z/V It, 


P = 


1 - 0)0 
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(6a) 

(6b) 









with the kinematic viscosity 


i/ = c2(r-l/2)5i, (7) 

where = c^/3. The above equation connects the fluid property to the 
parameter of LBGK model. 


3. incompressible D3Q14 MRT model 


Based on the above iD3Q15 LBGK model and the D3Q15 MRT model 
proposed by d’Humieres et al. [2^, we developed a three-dimensional 14- 
velocity incompressible (iD3Q14) MRT model. This model adopts the dis¬ 
crete velocity directions as 


{Cl, C2,. . . , C14} 


^ 1-10 0 0 0 1-1 1 -1 

001 - 10011 - 1-1 
^ 00001-111 1 1 


1 -1 1 -1 ^ 
1 1 - 1-1 
-1 -1 -1 -ly 


c, 


which do not contain the discrete velocity direction Cq, and we suppose c = 1 
such that the relevant quantities are dimensionless in the following. 

The evolution equation of iD3Q14 MRT model is 

14 

fa{x + Co,6t,t + St) - fa(x,t) = - ^ (/*(*, f) - t)), 0 = 1 - 14, 

i=l 

( 8 ) 

where is defined in Eq. ([2]) and Aq,* is the element of a 14 x 14 

collision matrix A. The above equation can also be written in a vector form: 


\f{x + cJt, t + St)) - \ f{x, t)) = -A{\f{x, t)) - t))), (9) 





where \f{x, t)) = {fi{x, t), f 2 {x, t), • • •, fuix, t))' is a column vector, and the 
superscript ' represents the transpose operator. For the iD3Q14 MRT model, 


we have defined a 

14 X 

14 transformation matrix: 







1 1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 \ 


-4 

-4 

-4 

-4 

-4 

-4 

3 

3 

3 

3 

3 

3 

3 
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1 

-1 

0 

0 

0 

0 

1 

-1 

1 

-1 

1 

-1 

1 

-1 


-4 

4 

0 

0 

0 

0 

1 

-1 

1 

-1 

1 

-1 

1 

-1 


0 

0 

1 

-1 

0 

0 

1 

1 

-1 

-1 

1 

1 

-1 

-1 


0 

0 

-4 

4 

0 

0 

1 

1 

-1 

-1 

1 

1 

-1 

-1 


0 

0 

0 

0 

1 

-1 

1 

1 

1 

1 

-1 

-1 

-1 

-1 

T = 
















0 

0 

0 

0 

-4 

4 

1 

1 

1 

1 

-1 

-1 

-1 

-1 


2 

2 

-1 

-1 

-1 

-1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

1 

1 

-1 

-1 

0 
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0 
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0 

0 

0 

0 

0 

1 

1 

-1 

-1 

-1 

-1 

1 

1 


0 

0 

0 

0 

0 

0 

1 

-1 

1 

-1 

-1 

1 

-1 

1 


1 0 

0 

0 

0 

0 

0 

1 

-1 

-1 

1 

-1 

1 

1 

-1 / 


( 10 ) 

to transform the distribution function into the moment with the linear map¬ 
ping m = T|/) and and simultaneously convert the colli¬ 

sion matrix into a diagonal one by A = TAT^^. Thus, we can further write 
Eq. ([9]) as 

\f{x + cJt, t + 6t)) - \f{x, t)) = -T~^A(m(a;, t) - t)), (11) 

where 

(^1 jX) Qx) jyt Qyi j Z) Qz) ^Pxx) Plouji Pxyi Pyzi Pxzi^xyz) ) (1^) 
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and the equilibria of the moments are 


p{eq) ^ Yp/S + U^/3, 
^{eq) = _7p ^2^ 


■(<^q) 

Jx = Mx, 
■(^q) 

Jy = Uy, 

■(eg) 

Jz = U,, 


Qx''^ = -7ux/3, 
= -7uy/3, 
= -7 uJ3, 


Wxx = 3ui - u% 

(eg) 2 2 

VuJU = %-u^^, 


(eg) _ 

Pxy — UxUy, 
(e?) 

Pyz = UyU^, 
(eij) 

Pxz' = UxU^, 

t^l = 0 . 


(13a) 

(13b) 

(13c) 

(13d) 

(13e) 

(13f) 


It should be noted that there are some differences in the dehnition of moments 
for iD3Q14 MRT model and d’Humieres et ah D3Q15 MRT model. Firstly, 
to construct a 14-velocity MRT model, we discard the moment related to 
kinetic energy square, which does not affect the hydrodynamics signihcantly, 
but is dehned in d’Humieres et ah D3Q15 MRT model. Secondly, based on 
the iD3Q15 LBGK model, we dehne a moment P, instead of the moment p. 


which is used in d’Humieres et al. MRT model. From Appendix A we will 
see that P is also a conserved quantity. 
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The diagonal collision matrix A is 


A = diag{si, S 2 , S 3 , S 4 , S 5 , Sg, S 7 , Sg, Sg, Sio, Sn, S 12 , S 13 , S 14 ), (14) 


where Si, S 3 , S 5 and S 7 are the relaxation parameters corresponding to the 
conserved moments. It should be noted that, the values of these parameters 
do not affect the recovered macroscopic N-S equations, which are set to be 
1.0 in the following simulations. In addition, Eq. (ITT)) can also be written as 


A - Sg, Sg, Sg, Sg, Sg, Sg, Sg, Sj/, S^, S^, S^y, Sj/, s^). 


(15) 


where Sg, Sg, Sg, s^, and s* are the parameters corresponding to the conserved 
moment, the moments related to kinetic energy, energy flux, viscous stress 
tensor and a third-order moment. This form of collision matrix is used in 
the Appendix A for recovering the incompressible N-S equations. Finally, it 
should be emphasized that iD3Q14 MRT model can recover to the iD3Q15 
LBGK model by setting A = (l/r)I, where r is the relaxation time of iD3Q15 
LBGK model and I is the unit matrix. 

The transformation matrix T can be obtained by two ways. The hrst one 
is similar to that of d’Humieres et ah MRT model, which uses the Gram- 
Schmidt orthogonalization procedure. Since this procedure for DdQ(g — 1) 
lattice model in three dimensions has not been given before, we show the 
details here. The transformation matrix T is obtained by orthogonalizing 
the matrix 

T = [CaxClyCazluxii, m,n,l> 0 , (16) 

where the element is a polynomial of m,n and I are integers. 

From a physical viewpoint, 14 row vectors of T correspond to 14 moments 
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of different orders. We have chosen T as 




( 17 ) 


where 


|0i)q: 11^0 11 ) 

I 4 ^ 2 ) Cl II II 1 

| 03 )a Cq-j;, 

|05)q ('ayi ^ 

107 ) 0 ; Cq, 2 , 

104)0 II^qII ^ox) 

|06)q 11^0 II 

|08)o 11^0 II ^oz) 


(18a) 


(18b) 


(18c) 


109)0 — 3c^^, 

I0io)o = cly - cL, 

1011) o CaxCayi 

1 012 ) 0 CayCdzi 

1013 ) 0 ('ax^-azi 

| 014 )q ^ax^ay^azi 

CK G {1, 2, • • •, 14}, ||Ca|| = {c^x + ^ay + Corresponding moments 

for 14 row vectors are the conserved moment (mi = P), the kinetic energy 
(m 2 = e), the momentum (m 3 , 5,7 = jx,y,z): ^^e energy flux (m 4 , 6,8 = Qx,y,z)^ 
the viscous stress tensor (mg = mio = = Pyy — mii,i 2 ,i 3 = 

Pxy,yz,xz) an third-order moment mi 4 = txyz- These moments can be 


(18d) 

(18e) 

(18f) 
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classified into four types according to the order. Among them, P is a zeroth- 
order moment, jx,y,z hrst-order ones, e and Pxx,ujLo,xy,yz,xz second-order 
ones, qx,y,z txyz are third-order ones. However, in the matrix T, 14 row 
vectors are organized in the order of corresponding tensors, rather than in the 
order of corresponding moments. The hrst two rows of T correspond to P and 
e, which are scalars or zeroth-order tensors. The next six rows correspond to 
Jxy Qx^ jyi %i 3z which are vectors or hrst-order tensors. The following 

hve rows represent the viscous stress, which is a second-order tensor. The 
last row represents a third-order tensor. After the orthogonalization of T, we 
obtain the transformation matrix T and the corresponding moments, which 

have the similar physical meanings. 

The above matrix T can be explicitly written as 

11111111111111^ 
1111 1 133333333 

1 -1 0 0 0 0 1 -1 1 -1 1 -1 1 -1 

1 -1 0 0 0 0 3 -3 3 -3 3 -3 3 -3 

0 0 1 -1 0 0 1 1 -1 -1 1 1 -1 -1 

0 0 1 -1 0 0 3 3 -3 -3 3 3 -3 -3 

0 0 0 0 1 -1 1 1 1 1 -1 -1 -1 -1 

, (19) 

0 0 0 0 1 -1 3 3 3 3 -3 -3 -3 -3 

33000033333333 
0011 -1 -1 00 0 0 0 0 0 0 
0 0 0 0 0 0 1 -1 -1 1 1 -1 -1 1 

0 0 0 0 0 0 1 1 -1 -1 -1 -1 1 1 

0 0 0 0 0 0 1 -1 1 -1 -1 1 -1 1 

0 0 0 0 0 0 1 -1 -1 1 -1 1 I -I J 

and T can be obtained by orthogonalizing the 14 row vectors of T in an 
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order. Supposing T = (|0i), |02), • • •, |014))^ then we have 


|02)q 


10l) Q II II ) 

7||c,||V2-15||c„||V2, 


|03)a (^axi 
|05)q Co-jy, ^ 

|07)« Caz) 

|04)q = (5 ||C«||^ - 13)Caa:/2, 
|06)a = (5 ||Ca||^ - 13)Cay/2, > 
|08)a = (Slicin' -13)c„,/2, ^ 

|09)a = 34 ^ - ||C„||% 1 

|0io)a = Cly - cL, j 


1011) o 
|012)o 

|013)a 

|014)o 


^ax^ayf ^ 
— ^ay^az’) 
^ax^az’) 

/ 


(20a) 

(20b) 

(20c) 

(20d) 

(20e) 

(20f) 


The first four orthogonal vectors are related to P and x-, y- and ; 2 -nionientnni 
modes: |0i)„ = 100^., |03)a = | 03 )«, I05)a = | 05 )q a^d | 07 )a = | 07 )c,. 
The vector |02)a is constrncted by orthogonalizing the energy mode |02)(j- 
Similarly, vectors |04 )q, |06)a and |08)a are respectively derived upon |04)^, 
|06)^ and |08)^. | 09 )q is built upon 109 )^ and |0io-i4)a = | 0 io-i 4 )q- ^ 
should be noted that, the row vectors in T are mntually orthogonal, bnt they 
are not normalized to assure that the elements of row vectors are integral, 
which can simplify the computation. Finally, from above derivation, we can 
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see that the orthogonalization procedure for DdQ{q — 1) lattice model is the 
same with that of DdQg lattice model, which promotes the generation of the 
second way to obtain the transformation matrix T. 

The second way to obtain the transformation matrix T is more straight¬ 
forward. This way properly uses the orthogonalization procedure, which has 
been done by d’Humieres et al.. The steps are as follows. First, we dis¬ 
card the hrst column and the third row of the transformation matrix for 
d’Humieres et al. D3Q15 MRT model, which means the distribution func¬ 
tion fo{x, t) in discrete velocity direction Cq and the moment e related to the 
kinetic energy square are not used. Secondly, to make the new 14 row vectors 
orthogonal to each other and assure the elements are integral, we can obtain 
the transformation matrix T in our model. 

ID3Q14 MRT model is described in details above. Then through the C-E 
expansion, we can prove that the above model can recover to the incom¬ 
pressible N-S equations as Eq. (E]) in the low Mach number limit with the 
kinematic viscosity 

= c2(r - l/2)5i, (21) 

where r = l/s^, = 1/3 (see Appendix A for details). 

4. Numerical results and discussion 

In this section, to verify the accuracy and stability of the proposed iD3Q14 
MRT model, the three-dimensional (3D) Poiseuille flow and pulsatile flow in 
a square channel and the 3D lid-driven cavity flow are simulated. In the 
simulations, the relaxation parameters in the collision matrix A are chosen 
to be Sc = 1.0, Se = 1.19, Sq = 1.2, St = 0.98, = 1/r, where r is related to 
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the viscosity by Eq. fl^TD . For the boundary condition imp lementation, the 


non-equilibrium extrapolation scheme is always applied 


26| . In addition, c is 


not always set to be 1 in the simulations. In this situation, u and p are hrstly 
normalized with c and and then used to calculate. After the computation 
is hnished, u and p are multiplied by c and as the hnal results. The 
relationship between the viscosity and relaxation time is the same with Eq. 
dH]), but Cg = c^/3. 


4-1. Steady flow: 3D Poiseuille flow 

The illustration of three-dimensional Poiseuille flow in a square channel is 
shown in Fig. [H The origin O is located at the center of the entrance plane. 
The flow region is dehned in a rectangular region: 0 < x < /, —a < y < a 
and —b<z<b, where I = 2,a = b = 0.5. Given the boundary condition: 


u(x, ±a, z, t) = u(x, I/, ±5, t) = 0, 


(22a) 


= Pin, p{l,y,z,t) = Pout, (22b) 


where Pin and Pout are the pressure at the entrance and exit of the channel. 


the three-dimensional Poiseuille flow has a steady solution [27|, 


u{y,z,t) 


16a2 

Z/TT^ 


dp, cosh{mz/2a), cos(i7iy/2a) 

/) 5 ^ (- 1 )^*')/'!- -3 ’ 

dx ^^ cosnltTTb/2a) 

(23a) 

v{x,y,zfl) =w{x,y,zfl) = Q, (23b) 


p{xfl) 


dp 

= Pin + 

dx 


(23c) 


where dp/dx = {pout —Pin)/l is the pressure gradient in the channel, v is the 
kinematic viscosity of fluid. 
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In the simulations, the initial and boundary conditions are set as 



( 24 ) 


where pin = 1.1 and pout = 1-0. To get the numerical solution, the criterion 
of steady state as follows is used. 


XI \u{xi,t + 5t) - u{xi,t)\ 


i 


< 1.0 X 10"^° 


(25) 




where X denotes the summation over the entire system. 


The simulations were carried out with the grid of 65 x 33 x 33 for u = 0.03. 
Fig. |2] shows the variation of u with y for different 2 ; location at section x = 1 
and the variation of p with x for different y location at section z = 0. To 
demonstrate the numerical results are independent of r, different r (1 /t = 
0.8,1.0,1.2) are used to simulate the flow. It can be seen from Fig. [2] that, 
the numerical results are in excellent agreement with the analytical solutions 
for different r. 

In Fig. [21 the velocity and pressure are taken from a: = 1 and z = 0 
sections, respectively. To validate the agreement is independent of sections, 
the variation of u with y at different x sections and the variation of p with x 
at different 2 ; sections are shown in Fig. [31 It can be seen that, the numerical 
results are still in excellent agreement with the analytical solutions. 

We also carried out the simulations with different lattice spacings. The 
global relative error in velocity held {GREu) between the numerical result 
and the analytical solution is dehned as 


GRE^ 



'.n - UaY + {Vn “ Va^ + (Wn “ Wa)^] 


(26) 
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where r;„, Wn, and Va-, Wa denote the numerical and analytical velocities, 
respectively, and the summation of i is over every grid point. The GREu 
for different lattice spacings at different 1/r are displayed in Table [H The 
variation of GREu with the lattice spacing is also shown in Fig. 01 The 
slopes of htting lines are about 1.85, 2.00 and 1.80 for 1/r =0.8, 1.0 and 1.3, 
respectively. This shows that the proposed iD3Q14 MRT model is of second 
order accuracy in space when simulating steady 3D Poiseuille flow. 


4-2. Unsteady flow: 3D pulsatile flow 

The 3D pulsatile flow is used to validate the proposed model for unsteady 
flow. The geometric conhguration of 3D pulsatile flow is the same with that 
of 3D Poiseuille flow, but the flow is driven by a periodic pressure gradient 
between the two ends of the channel. 


It should be noted that the 3D pulsatile flow in a circular pipe |28|, l29| is 


usually simulated by the lattice Boltzmann method, while the 3D pulsatile 
flow in a rectangular channel is little done before. We found an analytic 


solution 
paper 


or the 3D pulsatile flow in a rectangular channel in a very early 


301 1 . We hope this flow and its analytic solution can be widely used 


for validating the 3D LB models in the future. 

Supposing that the flow in the rectangular channel is laminar and incom¬ 
pressible, then the incompressible N-S equations for this flow are reduced 
to 


du 

'm 


dp d'^u 
dx^ ^ dy"^ 


d'^u 

dz^ 


1 


with the following pressure gradient imposed on the flow 


(27) 


^ = Gcos{ujt), 
dx 
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( 28 ) 











where G is an amplitude and a; is a frequency. The analytical solution of 


above flow is j30| : 


Re \ 

Uj 

f 1 O (-1)" rcosh(7„y/fe)cos(pnz/6) , ] 

1 ^ Pn I- cosh( 7 na/fo) 1 

n=0 

1 


a; 1 

cosh((T„^/b) cos(q„j//fe) 1 


\ 

1 1 

( cosh((T„) J J 

1 J 


(29) 


where 


7n = + (Tn = 


2 n + 1 

Pn = -Z-TT, 


2 ? 7 , + 1 b 

Qn = -^-TT-, 

2 a 


(30a) 

(30b) 


rj = b^^/ujju is the Womersley number. In the present study, the flow region 
is defined in a square channel with a = 6 = 0.5. 

The simulation parameters are set as follows. The grid size is Nx x 
Ny X Nz = 81 X 41 X 41, the period of the changing pressure is T = 100 
{u = 2 tt/T), and the magnitude of total pressure drop along the channel 
is Ap = 0.001 (G = Ap/l,l = 2), and the pressure at the outlet {pout) is 
set to be 1.0. St = 0.0125 is fixed in the simulations so that one period 
contains integral step. The initial state of velocity field is always set to be 
zero everywhere while the initial state of the pressure field is always set to 
be {pin+Pout)/ “2, where Pin = Pout — Ap is the pressure at the entrance (pj„ is 
determined by the Eq. (125]) at t = 0). The calculation of velocity field always 
began with several periods of initial steps to attain convergence criterion: 

X) \u{xi,t + T) - u{xi,t)\ 

(31) 


^ \u{xi,t + T)\ 


< 6.0 X 10 


-14 


where | • | denotes the absolute value operator and the summation of i is over 
the entire system. 
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We first carried out a set of simulations to get the velocity profiles across 
the channel at different times. Fig. |5] shows the velocity prohles in the lines 
X = 1, z = 0 at four different times: t = T/4,T/2,3T/4 and T. The relax¬ 
ation time r is chosen to be 0.6178 and 0.5500 with the Womersley number 
1 ] = 2.8285 and 4.3416, respectively. The velocity has been normalized with 
Umax, which is the maximal horizontal velocity of Poiseuille flow with pres¬ 
sure gradient —G. Umax = 1.876 x 10“^ and 4.420 x 10“^ for rj = 2.8285 and 
4.3416. The velocity profiles at different x and x locations are also plotted 
in Fig. inland Fig. [71 From above figures, it can be seen that the simulation 
results agree with the analytical solutions exactly. 

Next, we tried to get the accuracy of iD3Q14 MRT model for unsteady 
pulsatile flow by measuring the convergence order of GRE^- In the simu¬ 
lations, r and v are kept unchanged but c increases to two times when 5x 
decreases to half, rj = 2.8285. The simulation results at the time t = T/4, 
T/2, 3T/4 and T are shown in Table HI We also depicted the variation of 
GREu with the lattice spacing in Fig. [HI It can be seen that —ln{GREu) 
changes linearly with the —ln(Aa;). The slopes of fitting lines are 1.98, 1.87, 
1.98 and 1.86 for the solutions at the time T/4, T/2, 3T/4 and T, respec¬ 
tively. The above results show that the proposed iD3Q14 MRT model is of 
second order accuracy in space when simulating unsteady pulsatile flow. 

Finally, to compare the present MRT model with d’Humieres et al. MRT 
model, we calculate GREu at the time t = T for two MRT models. It should 
be noted that d’Humieres et al. MRT models have two types, one is based 
on Qian LBGK model and the other one is based on He-Luo LBGK model. 
Since He-Luo LBGK model has smaller compressible effect, we choose the 
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He-Luo type MRT model for comparison. In the comparison, the nnmerical 
results are obtained with the following parameters: 6t = 0.05 is hxed, and 
the grid is 41 x 21 x 21, c is equal to 1. r is set to be 0.9 while the relaxation 
parameter in d’Humieres et ah MRT model is set to be 1.0. Table [3]shows 
the variation of GRE^ with the pressure drop along the channel. It can be 
seen that the GREu of the present iD3Q14 MRT model is always smaller 
than that of d’Humieres et al. D3Q15 MRT model, which demonstrates 
that iD3Q14 MRT model is more accurate than d’Humieres et al. D3Q15 
MRT model for unsteady pulsatile flow. This may be attributed to that the 
compressible effect of our proposed model is smaller than that of d’Humieres 
et al. MRT model. 


4-3. 3D lid-driven cavity flow 

Because the 3D lid-driven cavity flow presents a variety of vortex motions, 
it is usually used in the validation of numerical method. Ku has simulated 
this flow with the pseudospectral method [3l|, and the result is widely used 
as a benchmark for the 3D lid-driven cavity flow. In this subsection, we also 
use the result by Ku to validate our method. 

The schematic of three-dimensional lid-driven cavity flow is shown in Fig. 
E We can see that the flow is conhned in a cubic box [0,1] x [0,1] x [0,1] 
and driven by the top lid, which is sliding with constant velocity Uo = 1.0. 
The flow in the cavity is supposed to be governed by the three-dimensional 
incompressible N-S equations. The Reynolds number is dehned as UoL/u, 
where L = 1.0 is the length of cubic box, u is the kinematic viscosity. 

In the simulations, the initial states of velocity and pressure helds, the 
velocities on the solid walls are all set to be zeros. At the section oi z = 0.5, 
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the symmetric boundary condition is set, i.e., 



( 32 ) 


At x = 0 and a; = 1 on ?/ = 1, u is set to be zero. At the points next to 
a; = 0 and a; = lon?/ = l,M = 0.3 and 1.0 are set. The setting of initial and 
boundary conditions is the same with that of Ku. 

In addition, c = 10 is hxed to make the flow satisfy the low Mach number 
limit. The convergence criterion for the velocity held is. 


XI X \uk{xi, t + 10005f) - Uk{xi, t) I 


i k 


< 1.0 X 10"^^ 


(33) 


XX + 1000(5t)| 


i k 


where the summations of i and k are carried out over all the grid points and 
the velocities in all directions. 

The grid independence of simulation result for iD3Q14 MRT model is 
hrstly examined. The Umin in the vertical center line [z = 0.5 and x = 0.5) 
for cavity how at Re=1000 by using four grids (49 x 49 x 25, 65 x 65 x 33, 
81 X 81 X 41 and 97 x 97 x 49) are -0.2619, -0.2693, -0.2730, -0.2751, 
respectively. The deviation of Umin between two adjacent grids are 0.0074, 
0.0037 and 0.0021. Taking Umin at grid 97 x 97 x 49 as a benchmark, the 
relative errors for three smaller grids from 49 x 49 x 25 to 81 x 81 x 41 are 
5.04%, 4.24% and 2.83%. 

Then the simulation results at grid 97 x 97 x 49 are used to compare with 
those of Ku for cavity how at Re=100, 400 and 1000. The comparison is 
shown in Fig. [TT] It can be seen that the agreement between our results and 
those of Ku are excellent. 
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Next, the stability of iD3Q14 MRT model is compared with that of 
iD3Q15 LBGK model. We hrstly use the iD3Q15 LBGK model to per¬ 
form the simulation of cavity flow at Re=1000 with the grid 97 x 97 x 49. 
The horizontal velocity prohle in the vertical center line is shown in Fig. [T2j 
It can be seen that the result of iD3Q15 LBGK model agree well with the 
results of iD3Q14 MRT model and Ku. 

Since many experimental and numerical studies on 3D lid-driven cavity 


flow 


32| are focused on the case of Re=3200, we also choose this case to 


test the stability of iD3Q14 MRT model. The simulation results of iD3Q14 
MRT model and iD3Q15 LBGK model with the grid 97 x 97 x 49 are shown 
in Fig. [131 In the hgure, the velocity vectors on the yz plane at x = 0.5 
and t = 50000(5f are plotted for two models. It should be noted that, in 
the plotting, the velocities on the left half plane are obtained from those on 
the right half one by the symmetry transformation, and only one grid point 
in every two is shown. From Fig. [121 it can be seen that both models can 
reproduce the Taylor-Gortler-Like (TGL) vortices (two pairs) at the bottom 
of the plane. It is also found that the flow above the TGL vortices (especially 
at the two top corners) are different for two models. This may be due to that 
the instantaneous flow tends to be affected by the model parameters for the 
unsteady cavity flow at Re=3200. 

Finally, we use iD3Q14 MRT model and iD3Q15 LBGK model to simulate 
the cavity flow at Re=3200 with the grid 49 x 49 x 25. The computation is 
stable for iD3Q14 MRT model, while the computation blows up for iD3Q15 
LBGK model. The result of iD3Q14 MRT model is shown in Fig. [HI The 
hgure shows the velocity vector plot of yz plane at x = 0.5 and t = 250005f. It 
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should be noted that St at the current grid is twice that at the grid 97x97x49. 
From Fig. [TU it can be seen that two pairs of TGL vortices can still be 
seen at the bottom of the plane. The above comparison demonstrates that 
iD3Q14 MRT model is more stable than iD3Q15 LBGK model. This may be 
attributed to that the iD3Q14 MRT model has more adjustable relaxation 
parameters. 

5. Conclusions 

We have proposed 14-velocity and 18-velocity multiple-relaxation-time 
lattice Boltzmann models for three-dimensional incompressible flows and re¬ 
covered the incompressible N-S equations through the G-E analysis. These 
two models only use 14 and 18 discrete velocities in velocity space and thus 
14 X 14 and 18 x 18 transformation matrixes in moment space, which can 
reduce the storage and computation costs in simulations. New models are 
constructed based on the three-dimensional He et ah LBGK models and 
d’Humieres et ah MRT LB models by realizing that He et al. LBGK models 
only need 14 and 18 discrete velocities without Cq = (0, 0), and at the same 
time, d’Humieres et al. MRT LB models contain a moment which has no 
effect on the recovered macroscopic N-S equations. 

We also have carried out the numerical simulations for the 3D steady 
Poiseuille flow, unsteady pulsatile flow in a square channel and the 3D lid- 
driven cavity flow using the proposed 14-velocity MRT model. For the 
Poiseuille flow, we have calculated the horizontal velocity prohles versus y 
at different x and z sections and the pressure prohles along the channel at 
different y and 2 ; sections. Our numerical results are in excellent agreement 
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with the analytical solutions. As to the pulsatile flow, the horizontal veloc¬ 
ity profiles versus y in the center line of the channel {x = 1 and z = 0) at 
four different times t = T/4, T/2, 3T/4 and T were computed. The velocity 
profiles at different x and sections were also measured. All these numerical 
results precisely agree with analytical solutions. For the lid-driven cavity 
flow, the velocity profiles in the vertical and horizontal center lines at section 
2 : = 0.5 were calculated for Re=100, 400 and 1000. The simulation results 
agree very well with the previous numerical results, which were obtained with 
an accurate pseudospectral method. 

For the steady Poiseuille flow and the unsteady pulsatile flow, we have also 
conducted the simulations to explore the numerical accuracy of the proposed 
MRT model. It is found that the proposed model has second-order accuracy 
in space. We also computed the global relative error in the velocity held of 
pulsatile flow, versus the pressure drop along the channel for the proposed 
MRT model and d’Humieres et ah MRT model. It is found that the global 
relative error of our model is always smaller than that of d’Humieres et ah. 
For the cavity flow, we have simulated the flow at Re=3200. It is observed 
that our MRT model can capture the TGL vortices at the grid 49 x 49 x 25, 
while He et ah LBGK model diverges at this grid. 

In short, we have developed two three-dimensional MRT LB models, 
which can recover to the incompressible N-S equations in the low Mach num¬ 
ber limit. These two models have higher storage and computation efficiency 
than the existing three-dimensional MRT LB models. The new models are 
based on He et ah LBGK models and d’Humieres et ah MRT models, but 
the new models are more stable or accurate. 
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Appendix A. The Chapman-Enskog analysis for iD3Q14 MRT LB 
model: derivation of incompressible Navier-Stokes 
eqnations 


In this section, we perform the C-E analysis for iD3Q14 MRT LB model 
to recover the incompressible N-S equations. With the knowledge that in the 
incompressible flow. 


0{6p) = 0{5p) = 0(M2), 

(A. la) 

0{u) = 0{M), 

(A.lb) 

where M represents Mach number, 6p and 6p are the pressure 

fluctuations, respectively. 

and density 

We hrst introduce the following expansions: 


OO 12 

fa{x + c„(5t, t = V —D'lfa{x, t), 

nl 

n=0 

(A.2a) 

OO 

n=0 

(A.2b) 

OO 

n=0 

(A.2c) 


where e = 6t and Dt = (d* + Cq, ■ V), we can rewrite the lattice Boltzmann 
equation 


14 

fc,{x + Co,St,t + 6t) - fa{x,t) = -'^A^i{fi{x,t) - f^^‘^\x,t)) (A.3) 

i=l 
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in different orders of e as follows: 


0(£“) : /f = /•'«>, 


(A.4a) 

(A.4b) 


i=l 


14 


14 


0(t") : ^ A„i/<‘>) = - 5 ; A„,/f >, (A.4c) 


where Dtp = + Cq, ■ V, and Eq. flA.4bp has been substituted into Eq. flA.4cP 

Then we can transform Eq. flA.4D into the moment space: 


^(0) _ 


(A.5a) 


(9tpl + Ckdk)m^°'' = 


(A.5b) 
(A.5c) 

Z, 

where = TC^T”^, Ck is a diagonal matrix with the k component of every 
discrete velocity Cq, as the element, and 


+ {dtol + Ckdk){l - = -kmP\ 


m 


A) = ('n p(") n n n nD) V -n = 1 9 

o , U, , U, (/y , U, (J^ , '^Pxx ’ ’ Pzx ’ ^xyz) ’ 


(A.6) 


It should be noted that, making the analysis to He et ah iD3Q15 LBGK 
model, which is similar to the analysis of Guo et ah LBGK model (see Eq. 

14 

(A15) in Ref. [l5|), we can deduce that rrii = Yl = (|p + + 

i=l 

0{e'^ + eM'^) = rri{^'’ + 0{e^ + eM‘^). This means that P is a conserved 
moment in the low Mach number limit, thus, P^^^ and P*^^^ are set to be 
zeros in Eq. flA.611 . 
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Expanding Eq. flA.Sbj) . we have 


7p/3 + u'^/Z 


Ux 


Uy 


Uz 


0 

— 7p + 


-5ua;/3 


— 5ii^ / 3 


—3uz(3 


see(l) 

Ux 


P + “x 


Ux Uy 


UxUz 


0 

-7uxt^ 


5ii^/3 - 411^ - 7p/3 


Ux Uy 


UxUz 


^qQx 

Uy 


Ux Uy 


P + ul 


UyUz 


0 

—7Uy/3 


UxUy 


5u^/3 - 4u^ - 7p/3 


UyUz 


^qHy 

Uz 

-\-dx 

UxUz 

-\-dy 

UyUz 

+8x 

p + 


0 




= — 

Sqq^z'^ 

-7ri^/3 


UxUz 

UyUz 


5u^/3 - ‘iul - 7p/3 


3u^ — 


Auxl3 


—2Uyj3 


—2uz/3 


3s^p£y 

2 2 


0 


2Uy j3 


—2uzl3 


Si/pL^2 

Ux Uy 


Uy /3 


Ux /3 


0 


(1) 

SuPxy 

UyUz 


0 


uz/3 


Uy /3 


SuPy]} 

UzUx 


Uzl3 


0 


Ux/3 


(1) 

SuPzX 

0 


UyUz 


UxUz 


Ux Uy 


1 ‘^tf'xyz 


(A.7) 


(A.7) 

From above eqnations, we can obtain 

dto (7p/3 + mVS) + dxUx + dyUy + = 0, (A.8) 

^to^x ^xi^P ^yip^x'^y) A d^(^UxUz) 0? (A.9a) 

^to^y A dxiUxUy) A ^y(^P A A dziUyUz) 0, (A.9b) 

"h ^xip^x'^z) A A dzij^ A *^ 2 ) (A.9c) 


Since 0{5p) = 0{M‘^) and 0{u) = 0{M), we have 5tg(7p/3 + u‘^/3) = 
0{M‘^). Omitting the 0{M‘^) term, Eq. flA.Sp becomes 

dxUx + dyUy + = 0, (A.IO) 

which is the continnity eqnation of incompressible N-S eqnations. 

From the expansion of Eq. flA.Scp . we can have 

dt^Ux + ^,,[(1 - Sy/2)p^xl + 2(1 - Se/2)e(^V21] + dy[{l - Sj./2)pS] 

+ 9^[(1 - Sy/2)p^zl] = 0 , 

(A.lla) 
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dt^Uy + 5a;[(l - Sy/2)p^xy] + dy[{l - sy/2)- p^xl)/2 + 2(l - Se/2)e^^y2l] 

+a,[(l-s,/2)pg] =0, 


(A.llb) 

dt,u, + 6>,[(1 - sj2)p^^I] + a,[2(l - Se/2)eW/21 - (1 - sj2){p^^l+p^^})/2] 

+a,[(l-s,/2)pg] = 0, 


Combining Eq. flA.91) and Eq. 
we can obtain 


(A.llc) 

flA.lip with the operation flA.9D +£X flA.lip . 


A dxip A A dy{uxUy^ A d^ipix^z^ / ^^P^y\ 

-edx[{l - s./2)pW + 2(1 - Se/2)eW/21] - £^,[(1 - s./2)pW], 


(A.12a) 


A dxi^UxIly^ Oy(^p ^^^y) ] 

-edy[{l - Sy/2){p\^l -pxI)/2 + 2(1 - Se/2)eW/21] - - Sy/2)p^y)], 

(A.12b) 

dtu^ + dx{uxu^) + dy{uyu^) + d^{p + ul) = -£^[(1 - s^/2)p^^J] 

-edy[{l - Sy/2)p^y}] - £^^[2(1 - Se/2)eW/21 - (1 - Sy/2){j)'^l +Pxl)/‘2]. 




(A.12c) 

According to Eq. flA.7P. yields 


e(i) = 

1 5 

[9^0 ( Tp -\- U ^ — (dx^x “1“ dyUy -\- ^2:^2:)]: 

Sg o 

(A.13a) 

= - 

fxx 

1 4 2 2 

2^ It ) + ^dxUx ^dzUz], 

(A.13b) 


- ul) + ‘^dyUy - ‘^dzUzi 

(A.13c) 


1 11 

Pxy ~ [dto{UxUy) + —dxUy + —dyllx^, 

Si/ o o 

(A.13d) 


p\iz ~\dtQ{UyUz) + —dyllz + ’^dz'Uy], 

(A.13e) 
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(A.13f) 


pi:x [^io(^2^x) “1“ cjdxUz “1“ Q^z^a:]) 

5;y O O 


Using Eq. flA.9p . we can estimate the terms dt^ul, dt^ui and dt^ul as 


^to^x 2'Ua;(5x(P A '^x) A dy(UxUy) dziUxUzy)-, 


(A. 14a) 


^Uyidxiux^'y) A dy{j)uy) ~|“ dziviyUz)')^ (A.14b) 

^to^z ‘^'^z{dx{UxUz') A dy{UyUz) A 9z (^P A ^ 2 :))' (A.14c) 

From Eq. flA.ip . we can see that dtoU^, dt^Uy and are all in the order 


of 0{M^). Omitting the 0{M^) terms, Eq. flA.13j) becomes 


e(i) = 


3 s, 


p(i) = 
Irxx 


95, 


'{dx'Ux 9? 


'{^dx'Ux dy^y 5 


p(l) = 
Iruju) 


pd) = 
t'xy 


pd) — 
fyz 


p(d = 

f zx 


3s, 


-(dylly dzUz')i 


3s, 


'(^dxUy A dyUx)^ 


3s, 


3s, 


-{dyllz A dzUy')i 

'{dxUz A dz,Ux)^ 


(A.15a) 

(A.15b) 

(A.15c) 

(A.15d) 

(A.15e) 

(A.15f) 


where Eq. flA.Sp has been substituted into Eq. (IA.15aP and the term dt^v? 
has been omitted. 

Substituting Eq. (IA.15P into Eq. flA.12p . and supposing z/ = c?fr—l/2l^t. 
where r = l/s„ and = 1/3, then we have 


dtUx A dxiul) + dyiuxUy) + d^iuxU^) = -dxP A v{dlux A d^Ux A d^Ux), 

(A.16a) 


30 
























dtUy + dx{Ua;Uy) + dy{ul) + d^{UyUz) = -OyP + + OyUy + d^Uy) , 

(A. 16b) 

dtUz + d^{u^u^) + dy{uyu^) + d^{ul) = -d^p + u{dlu^ + d'^u^ + 

(A.16c) 

which are the momentum equations of incompressible N-S equations. It 
should be noted that the continuity equation and e = 6t have been used in 
deriving the above equations. 

All in all, from the Chapman-Enskog analysis for iD3Q14 MRT LB model, 
we can see that, in the low Mach number limit, the proposed 14-velocity 
model can recover to the incompressible N-S equations, which can be written 
in a vector form as Eq. l|6|). 

Appendix B. iD3Q18 MRT LB model 

The three-dimensional incompressible MRT LB model with 18-velocity 
adopts the following discrete velocity directions: 

{Ci, C2, . . . , Cis} = 

1 -1 00001 -1 1 -1 1 -1 1 -1 00 0 0 

<001 -1 0011 -1 -1 00 0 01 -1 1 -1 

00001 -1 00 0 011 -1 -1 11 -1 -1 

\ 

Ui_Q = 1/18, a; 7 _i 8 = 1/36; = c^/3, c = 1. 

(B.l) 
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18 orthogonal basis vectors can be derived by the Gram-Schmidt orthogo- 
nalization procedure: 


101)« II II ) 

|02)a = 3 ||C„||^ - 5 ||C„||° 


l^s)^ Ca.xi 

1 ^ 5 ) 0 : ^ 

10?) a ^a.z^ 

/ 

1^4)0: (5 ||cq.|| 9)cq.3;, 

|06)a (5 ll^all 

\4^s)<y (5 ll^all 9 )Cq-^, 

1 ^ 9 ) 0 : = 3Cq,^ — ||cq,|| , 

\<Pn)a = cly-cl^, 


> 


\<Pl3)a 

\<Pl4)a 

|015)a 


^ax^ay") 
^ay^az") 
— ^ax^az-) 


> 


(B.2a) 


(B.2b) 


(B.2c) 


(B.2d) 


(B.2e) 


|0io)a = (3 ||c«||'-5)(3c: 
|012)a = (3 ||Ca||^ - 5)(c; 

|016)a = {cly - 
|017)o = (c; 


2 

ax 


2 

ay 


I Co 
^2 


-,2 

^az 


^axJ^ciyi 


-,2 

^ax 


^ay)^0!Z: 


(B.2f) 


(B.2g) 


I018)a = (Co 

where a G {1, 2, • • •, 18}. The corresponding 18 moments {m^(a;,t)l0 = 
1, 2, • • •, 18} are dehned as: 


Tn(^Xj t) (T*, e, jxj Qxi iy) Qyi iz) Qzi 3pxx) “i^^xx) Ploloi Pxyj Pyz) Pxzi tx: ty: 


tzY- 


(B.3) 
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The diagonal collision matrix is 


A. — d,icig(^S(;, Sg, Sg, Sg, Sg, Sg, Sg, Sg, s^,, Stt, Sj/, Stt, Sj/, Siy, Si,, Sj, Sj, Sj), (B.4) 
and the transformation matrix T can be obtained from Eq. (IB.2p . 


/ 1 1 1 1 

-2 -2 -2 -2 


T = 
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0 

0 

0 
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-1 

4 
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0 

0 

0 

2 


-4 -4 

0 0 


0 

0 

1 

-4 

0 

0 

-1 

2 

1 


0 

0 

-1 

4 

0 

0 


1 1 

-2 -2 

0 0 


0 

0 

0 

1 

-4 


-1 -1 

2 2 


1 


-2 -2 

0 0 


-1 

2 

0 

0 

0 

0 

0 

0 


0 

0 

0 

-1 

4 

-1 

2 

-1 

2 

0 

0 

0 

0 

0 

0 


1 

1 

-1 

-1 

1 

1 

0 

0 

1 

1 

1 

1 


1 

1 

1 

1 

-1 


1 

1 

-1 

-1 

-1 


-1 -1 

0 0 


-1 -1 

0 0 


0 

-1 


-1 -1 

0 0 


0 

1 

1 

1 

1 

1 

0 

0 

-1 

1 

0 


1 

1 

-1 

-1 
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0 

1 

1 

1 

1 


-1 -1 

-1 -1 

0 0 


0 

1 

-1 

0 
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The equilibria of 18 moments are defined as 

p(e9) = 2p + 

gleg) = _p _|_ 


0 

-1 

1 

0 
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1 

1 

1 

1 

0 

0 

-1 

-1 

1 

1 

-1 

-1 

0 

0 

-1 

-1 

0 

-1 


1 

1 

-1 

-1 

0 

0 

-1 

-1 

1 

1 

-1 

-1 

0 

0 

1 

1 
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1 

1 

0 

0 

-1 

-1 

1 

1 


1 \ 
1 
0 
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-2 -2 

-2 -2 

0 0 


1 -1 

1 -1 

-1 -1 

-1 -1 

-2 -2 

-2 -2 

0 0 


-1 -1 


0 

0 

-1 

0 

0 

-1 

-1 


-1 1 

0 0 

0 0 

1 -1 
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(B.5) 


1 / 


(B.6a) 


■(eg) 

Jx = Mx, 
■(eg) 

Jy = Uy, 

■(eg) 

Jz = U,, 
(eg) 2 

Qx = -gMx, 

(eg) 2 

q^y = --Uy, 

(eg) 2 

qz = 


(B.6b) 


(B.6c) 
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(B.6d) 


Spxx = Sul ~ 

= ul-ul, 


StTxx'^ = -liSul - u^) 


= -\{ul-ul), 


(eq) _ 

Pxy UxUy^ 

Pyz = UyU^, 
(eg) 

Pxz = UxU^, 


tf' = 0, 

4“" = 0, 

4'” = 0. 


(B.6e) 


(B.6f) 


(B.6g) 


Through the C-E expansion, the incompressible N-S equations can be 
recovered in the low Mach number limit for iD3Q18 MRT model. The re¬ 
lationship between the kinematic viscosity and the relaxation parameter Siy 
is the same with that of iD3Q14 MRT model. ID3Q18 MRT model can re¬ 
cover to the corresponding He et al. iD3Q19 LBGK model by setting all the 
relaxation parameters to be 1/r, where r is the relaxation time of iD3Q19 
LBGK model. 

It should be noted that, we also simulated the steady Poiseuille flow, 
unsteady pulsatile flow and lid-driven cavity flow in three dimensions with 
iD3Q18 MRT model. The numerical results show that iD3Q18 MRT model is 
also of second order accuracy in space for steady Poiseuille flow and unsteady 
pulsatile flow, but iD3Q18 MRT model is more stable than iD3Q14 MRT 
model for the lid-driven cavity flow. In the simulation of lid-driven cavity 
flow with the grid 49 x 49 x 25 , the reached Reynolds number for iD3Q18 
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MRT model is larger than that for iD3Q14 MRT model. This hnding is 


consistent with the observation that D3Q19 MRT model is more stable than 


D3Q15 MRT model by d’Humieres et ah 


24| . In addition, based on our 


construction method of {q — 1) x [q — 1) transformation matrix in three 
dimensions, we can also obtain iD3Q26 MRT model for three-dimensional 
incompressible flows. A detailed comparative study of above three MRT 
models is left for future work. 
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Figured) The schematic of three-dimensional Poiseuille flow. 

Figure [2) The left hgures show the horizontal velocity prohles for Poiseuille 
flow at section x = 1, while the right hgures show the pressure prohles at 
section z = 0. lines: analytical solutions, symbols: numerical results. 

Figure O The velocity and pressure prohles at diherent x and sections 
for Poiseuille how. 

Figure d) The variation of GREu with the lattice spacing for Poiseuille 
how by using diherent 1/r. 

Figure [S] The variation of u with y for pulsatile how at the location 
X = 1, 2 ; = 0. lines: analytical solutions, symbols: numerical solutions. 

Figure [6l The variation of u with y for pulsatile how at diherent z loca¬ 
tions at section x = 1 for r] = 2.8285. 

Figured) The variation of u with y for pulsatile how at diherent x loca¬ 
tions at section z = 0 for rj = 2.8285. 

Figure [S) The variation of GREu with lattice spacing for pulsatile how 
at diherent times, r] = 2.8285. 

Figure [9) The schematic of three-dimensional lid-driven cavity how. 

Figure dQ) The distribution of u at four diherent grids in the vertical 
center line (z = 0.5 and x = 0.5) for cavity how at Re=1000, (b) is the 
magnihcation of square area in (a). 

Figured!) The velocity distribution in the vertical and horizontal center 
lines at section z = 0.5 for cavity how at diherent Re, □ : the results of Ku, 
- : present simulation. 

Figured!) The distribution of u in vertical center line for cavity how at 
Re=1000. The simulations are carried out with grid 97 x 97 x 49. 
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Figure US) The velocity vector plot of yz plane aX x = 0.5 and t = 50000(5t 
for cavity flow at Re=3200, the grid is 97 x 97 x 49. The length of the arrows 
is three times the actual velocity magnitude. 

Figure fm The velocity vector on the yz plane at x = 0.5 and t = 250005f 
for cavity flow at Re=3200, the grid is 49 x 49 x 25. The length of the arrows 
is three times the actual velocity magnitude. The result is obtained from 
iD3Q14 MRT model, while the simulation by iD3Q15 LBGK model diverges. 


41 


Table 1: The GRE^ of Poiseuille flow for different lattice spacings and different 1/r. 




GREu 


1/r 

Ar = 1/8 

Ax = 1/16 

Ax = 1/32 

Ax = 1/64 

0.8 

1.116 X 10"^ 

3.277 X 10-2 

9.001 X 10-3 

2.371 X 10-3 

1.0 

2.976 X 10-2 

7.400 X 10-3 

1.846 X 10-3 

4.610 X 10-^ 

1.3 

5.824 X 10-2 

1.854 X 10-2 

5.232 X 10-3 

1.392 X 10-3 

ble 2: The GREu of pulsatile flow for different lattice spacings and different tin 



GREu 


t 

Ax = 1/20 

Ax = 1/40 

Ax = 1/60 

Ax = 1/80 

T/A 

1.750 X 10-2 

4.465 X 10-3 

1.994 X 10-3 

1.124 X 10-3 

T/2 

4.060 X 10-2 

1.169 X 10-2 

5.444 X 10-3 

3.134 X 10-3 

3T/4 

2.201 X 10-2 

5.661 X 10-3 

2.535 X 10-3 

1.430 X 10-3 

T 

3.811 X 10-2 

1.108 X 10-2 

5.173 X 10-3 

2.981 X 10-3 


Table 3: The variation of GREu with the pressure drop for pulsatile flow by using different 
MRT models, t = T. 


Ap 

Re 

Umax 


GREu 

d’Humieres et al. D3Q15 MRT 

iD3Q14 MRT 

0.001 

0.41 

0.0055 

0.0096 

0.0067 

0.0059 

0.005 

2.1 

0.028 

0.047 

0.0078 

0.0067 

0.01 

4.1 

0.055 

0.096 

0.0092 

0.0075 

0.02 

8.3 

0.11 

0.19 

0.012 

0.0090 

0.05 

20.7 

0.28 

0.48 

0.020 

0.010 
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Figure 1: The schematic of three-dimensional Poiseuille flow. 
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(a) horizontal velocity profiles at different x sections 


section: z=-0.125; 1/r=1.2 



y 

(b) pressure profiles at different sections 


section: y=-0.125; 1/r=1.2 



Figure 3: The velocity and pressure profiles at different x and z sections for Poiseuille 
flow. 
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Figure 6: The variation of u with y for pulsatile flow at different 2 ; locations at section 
x = 1 for 7? = 2.8285. 
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Figure 9: The schematic of three-dimensional lid-driven cavity flow. 
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- 0.32 - 0.3 - 0.28 - 0.26 - 0.24 - 0.22 - 0.2 - 0.18 - 0.16 - 0.14 


u 

Figure 10: The distribution of u at four different grids in the vertical center line (z = 
and X = 0.5) for cavity flow at Re=1000, (b) is the magnification of square area in (a) 
















Figure 12: The distribution of u in vertical center line for cavity flow at Re=1000. The 
simulations are carried out with grid 97 x 97 x 49. 
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(a) iD3Q14 MRT model 



(b) iD3Q15 LBGK model 



Figure 13: The velocity vector plot of yz plane at x = 0.5 and t = 50000<5t for cavity flow 
at Re=3200, the grid is 97 x 97 x 49. The length of the arrows is three times the actual 
velocity magnitude. 
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iD3Q14 MRT model 
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Figure 14: The velocity vector on the yz plane at a; = 0.5 and t = 25000(5t for cavity flow 
at Re=3200, the grid is 49 x 49 x 25. The length of the arrows is three times the actual 
velocity magnitude. The result is obtained from iD3Q14 MRT model, while the simulation 
by iD3Q15 LBGK model diverges. 
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